Freezing effects in the two dimensional one-component plasma and in thin film type II 

superconductors 
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We present results of Monte Carlo simulations of the two dimensional one-component plasma 
and of the Ginzberg-Landau model in the lowest Landau level approximation, with both charges 
and vortices respectively confined within a disc. In both models we see that as the temperature is 
reduced, oscillations in the radial density develop which spread into the bulk from the edge of the 
disc. The amplitude of these oscillations grows as the temperature is lowered and the length scale 
over which the oscillations occurs is the same as the correlation length for local crystalline order 
at that temperature. At temperatures similar to those where earlier studies have reported a first- 
order fiuid-crystal phase transition, the correlation length is comparable to the linear dimensions 
of the samples studied, which suggests that finite size effects will be affecting the accuracy of their 
conclusions. 



PACS numbers: 05.20.-y, 61.20.Ja, 74.62.Yb 

In this paper, we report on freezing-related effects in 
two models with long-ranged interactions between 'par- 
ticles', both of which have some bearing on the behav- 
ior of clean thin-film superconductors. The models are, 
the one-component plasma in two dimensions and the 
Ginzberg-Landau model in the lowest Landau level ap- 
proximation; henceforth 2D OCP and the GL LLL mod- 
els respectively. 

Published simulations of these models have given con- 
flicting results regarding even the existence of the fluid- 
crystal transition. For the 2D OCP with periodic bound- 
ary conditions, Leeuw and Perraro^ found a first order 
transition at about F = 140 (with that parameter defined 
below) which showed up as hysteresis in the energy when 
the system is cooled and then heated through the transi- 
tion temperature. The same conclusion was reached by 
Choquard and Clerouin" in simulations of the disc geom- 
etry with some particles close to the boundary fixed at 
vertices of a triangular lattice. In contrast, Monte Carlo 
simulations for charges on a spherical surface by Moore 
and Perez-Garridc"^ found no evidence for a transition; 
instead the correlation length for short-range crystalline 
order in the fiuid was found to grow to infinity as the tem- 
perature was reduced to zero. A dynamical investigation 
of the disc geometry, this time without any constraint 
on the particle positions has been made by Ying-Ju Lai 
and Lin I^. They looked at the dislocation density in dif- 
ferent regions and suggested that the system melts from 
the outside as dislocations proliferate around the intrin- 
sic disclinations which are close to the boundary at low 
temperatures. They found that the central region loses 
translational and then orientational order as the temper- 
ature increases. 

For the GL LLL model, Monte Carlo simulations for 
two-dimensional periodic cells^i^ii have supplied evidence 
for a first-order fluid-crystal transition at an effective 
"temperature" q;t = — 10 (ay is defined below). Sim- 
ulations on a spherical surface^ and within a diso^ find 
only a fiuid phase with a growing correlation length for 



decreasing effective temperature. 

Because whether one sees a phase transition or not 
seems to be dependent upon the boundary conditions 
used, we have revisited the problem to investigate how 
finite size effects may be affecting the outcomes of these 
diverse simulations. In particular, we have studied the 
radial density for the 2D OCP and the GL LLL mod- 
els when the particles are confined to a disc. We have 
found that it develops oscillations or ripples of increasing 
amplitude as the temperature is reduced which reach in 
from the edge of the disc by a distance of the order of the 
length scale of short-range crystalline order in the system 
at that temperature, as found for example in the work 
of Moore and Perez-Garrido2.. The existence of such os- 
cillations has been found in other fluid systems in the 
presence of an interfac o^°'^^ . 

Perhaps the most straightforward way of seeing the 
origin of the oscillations induced in the density by the 
presence of a surface is via linear response theory. A 
small external potential 5(j){r), whose Fourier transform 
is 50(k), will induce a change in the one-particle density 
Sp{r). Its Fourier transform is given by^^^ 



(5p(k) = -/3pS'(k)J0(k) 



(1) 



so the density changes Sp{r) induced by the potential will 
depend on the analytic structure in the complex k plane 
of the structure factor S(k) = (/3(k)p(— k)) /N of the bulk 
system. It is argued that at distances well away from the 
surface its effects can be treated as if it corresponded to 
a perturbation applied at the edge of the system. When 
there is local crystalline order, the structure factor has a 
peak close to the magnitude of the first reciprocal lattice 
vector, say at |Ko| and the profile of the peak is roughly 
Lorentzian 



S{k) 



1 



(fc~ifo 



1/e 



as seen in earlier work^,, with ^ the bulk correlation 
length, which has poles at k = Kq ± i/^. The imaginary 
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part controls the length scale of decay of the density os- 
cillations. At large distances d from the edge of the disc 



Sp{d) ~ exp(— d/^) cos{Kod). 



(2) 



We first introduce the models, and then present our 
results on the surface induced effects in a disc geometry. 
Finally we mention the implications of these results for 
numerical simulations of the 2D OCP and the GL LLL 
model around the reported phase transition. 

The 2D OCP is a system of particles interacting loga- 
rithmically in a central harmonic confining potential. We 
arrive at the harmonic confining potential naturally by 
finding the energy of interaction of each charge with a 
background charge density which is spread out over the 
whole plane. Suppose we have N identical charges e on 
a fixed background of uniform charge density which con- 
fines them within a disc of radius R. It can be shown that 
the energy V oi a configuration of N charges each labeled 
by integer i and with position {vi} is, up to constants. 
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(3) 



The temperature parameter that we use is F = /kT. It 
is convenient to set the average number density N/ ttR^ to 
I/tt. In the London regime of thin film type II supercon- 
ductors, which is characterized by low fields and hence 
low vortex densities so that variations in the complex 
order parameter can be neglected, vortices also interact 
logarithmically as in the 2D OCP. 

The Ginzberg-Landau free energy functional for su- 
perconductors at temperature T depends on a spatially 
varying, complex order parameter ip (r), which is coupled 
to a magnetic vector potential A. The superconductor is 
exposed to a uniform and static magnetic field B in the 
C direction. This field is close to the upper critical field 
(deduced from the mean field approximation) so that it 
is almost uniform within the sample; we neglect the con- 
stant free magnetic field energy term. 



F= f d^rl—\(-ih\/ -e*A)i;f 
J 1 2m 



The integral is taken over the entire sample volume. The 
parameters a, f3, and m are related to the temperature 
dependent penetration depth and correlation length of 
the superconductor. We shall be interested in a sample 
that is sufficiently thin (thickness t) in the direction of 
the field that we can neglect any variation in the order 
parameter in this direction and remove the integral over 
Then we restrict the order parameter to a linear com- 
bination of Landau levels of lowest energy. In the sym- 
metric gauge, with complex coordinates (A^ = —By/2 
and Ay — Bx/2), these take the form 



where I = 0,1,... and gi = (Trm!)"^/^ (1/2P)('"+^)/^ 
are normalization coefficients. P = h/ B is the squared 
magnetic length. The order parameter is 



N 

'E 



vi4>i (z) ■ 



where Q = (27rP//?t)^/^ is introduced now to obtain a 
neater expression later. The sum runs over the iV -I- 1 
lowest angular momentum states. This is a polynomial 
of degree A^ in z which therefore has A^ zeroes correspond- 
ing to the vortex positions in our finite system. This is 
substituted into the free energy functional. When the 
integrals have been performed we get 
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0, (z) = 5^^'e-N^/4^, 



p,q,r,s=Q V-F-y- ■ -J 

where ar — Q^t{a+ {e* Bh/m))^^. Within this approx- 
imation, the vortex matter has identical properties along 
each of a family of curves labeled by ax in the magnetic 
field - temperature plane. The quantity exp {—F) is taken 
as a Boltzmann weight in the calculation of the canoni- 
cal partition function Z. Several authors^"'' have shown 
that at sufficiently high field in Type II superconductors, 
ffuctuations in the magnetic field can be neglected. 

As for our simulations, we used the Metropolis 
algorithmic to approximate thermal averages by mak- 
ing pseudorandom steps in the coefficients {vi} in the 
Ginzberg-Landau simulations or in the positions of the 
charges in the OCP simulations. We thermalized the 
sample for a period and then started taking measure- 
ments of configurations to contribute to the averages. 

We computed the radial density for the 2D OCP as 
the thermal average of a fine-grained histogram of parti- 
cle numbers in circular shells for different discrete values 
of r. As is known from the exact solution for F = 2, the 
radial density is uniform from the center to the edge at a 
value of I/tt. At lower temperatures, oscillating behavior 
creeps in from the edge of the sample. With A^ — 400 the 
pattern spreads across the system for F = 100 though the 
amplitude of the peaks continues to rise as the tempera- 
ture is lowered (Fig. [1]). We find, on increasing the num- 
ber of steps in the simulation that these features do not 
change, indicating that our system has been adequately 
equilibrated for their study. 

We first consider the outermost peak. For F greater 
than about 100 this peak corresponds well to an isolated 
shell (with a fixed number of particles). This is because 
the density on either side of the outer peak drops almost 
to zero, so there must be very little particle exchange 
between this peak and neighboring ones. We have found 
that the height of this peak is proportional to vT. The 
reason for this is as follows. The radial oscillations of 
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FIG. 1: Radial particle densities for the 2D one-component 
plasma at T = 20, 40, 60, 80, 100, 120 and 140 for 400 charges 
showing progressive ordering from the edges. We have sepa- 
rated the curves along the vertical axis for the sake of clarity. 



particles in this shell have an amplitude X that is equal 
to half the width of the outer shell. The average energy 
of this mode e^(X^) is, by equipartition, of order kT. So 
the mean displacement of the particles from equilibrium 
is roughly proportional to I/Vt and hence, the height 
of the peak is proportional to Vt. The other peaks do 
not obey this relation as closely as the outermost peak 
because the rest of the sample has roughly triangular 
order and so the assumption of isolated shells with radial 
oscillations breaks down. 

It is natural to expect that the peak heights in the 
radial density will decay with distance from the edge of 
the disc on a length scale which is a reflection of the 
short-range crystalline order in the system and that the 
length scale involved is just the correlation length f of 
the 2D OCP. Moore and Perez-Garridc'^ have studied 
the structure factor S'(k) of the 2D OCP on the surface 
of a sphere and found that from the inverse width of its 
first peak that the correlation length is ^ y/V at low 
temperatures: that is, when F > 60. We make the as- 
sumption that the same behavior for occurs in our disc 
geometry, h — Iiq is the radial density peak height above 
the uniform radial density ho which is ~ I/tt. We plot 
ln[(/i — /io)/Vr] against r/^ where r is the distance from 
the outer peak where the ^/T scaling is based on the vari- 
ation of the absolute height of the outermost peak. The 
correlation length is taken to vary as VT — P where (3 is 
a fitting parameter. Fig. [5] shows collapse of data at var- 
ious temperatures onto a straight line revealing that the 
peak heights vary exponentially with distance from the 
highest peak with a characteristic length that varies with 
temperature in the same way as on a spherical surface for 
F > 60. The parameter /3 that gives the best collapse is 
also consistent with already published data for charges 
on a spherical surface^. 

Because the rather general considerations of Evans and 
coworkers give oscillating exponential behavior as in 
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FIG. 2; Plot oi {h ~ ho)/VT against r/^ where h are peak 
heights in the radial density, a is the mean value of the ra- 
dial density and r is the radial distance from the edge of the 
disc (as measured from the maximum peak). The correlation 
length ^ = \/r — (3 with fitting parameter f3 — 4.46. Plot 
includes data for F = 80, 100, 120 and 140. 
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FIG. 3: Two plots of radial densities for different system sizes 
(100, 400, and 600 charges) superimposed with the edge peak 
superimposed. The plot shown is for F = 100. 



Eq. [5] only asymptotically at large d, we should not be 
surprised that at least the first peak height is exceptional 
in Fig. [2j nor that the peaks close to the edge have 
varying width and shape. Indeed, if this relation were to 
hold close to the edges with \/T scaling of the heights, the 
density would not be completely positive even at quite 
high temperatures. We expect the oscillating exponential 
to be a good fit to our data only asymptotically. 

The data in Fig. [2] were taken for a disc containing 
400 charges. We have superimposed 2D OCP radial den- 
sities for identical temperatures and different system sizes 
with the outer peak as a reference point (Fig. [3|). The 
superimposed plots overlap almost perfectly independent 
of the system size, indicating that the oscillations them- 
selves are not finite size effects, but are instead related 
to the presence of an edge in the system. 

We now consider the thermal averaged superconduct- 
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FIG. 4: Order parameter radial densities for equal to —2, 
—4, —6, —8 and —10 from bottom to top. There are 400 
vortices in the disc. 

ing order parameter density {\ip{r)\'^) as a function of the 
distance from the center of the sample r in the GL LLL 
model. This can be computed from the N + 1 averages 
{vjnV^) as follows 

n=0 

When ax is positive, the order parameter density is uni- 
form except for a sharp drop to zero at the edges of the 
sample. Plots of the radial density for various negative 
values of ar and N — 400 are shown in Fig. |4l As 
Q!T is lowered, peaks form at the edges and increase in 
amplitude; shell-like structure spreads into the sample 
from the edges. At about ar = —8, we observe that the 
ripples have spread all the way across the sample. If we 
look at the particle configurations for negative ar, we see 
that the peaks are due to increasingly shell-like ordering 
towards the edges. We have confirmed (for small sys- 
tem sizes for which we can compute the vortex positions 



quickly) that the troughs in the radial density match up 
with regions of high vortex density as we would expect. 
The tendency to form a triangular crystal at low tem- 
peratures emerges quite close to the edge of the sample 
so that the evolution from shell-like order to triangular 
order happens within a few lattice spacings. There is no 
evidence for a phase transition to a crystal in this sys- 
tem; instead the ordering is gradual^. The period of the 
radial density ripples is the same for all system sizes and 
temperatures. In fact, if we compare the radial densities 
for samples of different TV at a given temperature, we 
find that they very nearly overlap close to the edges for 
N = 400 and N = 600, but the N = 100 discs show pro- 
nounced differences compared to the other two sizes. We 
have repeated the scaling analysis that we applied to the 
radial density peaks of the 2D OCP, this time with the 
assumption that ^ ^ (|q;t| — /3) as suggested for ar < — 8 
in work by Dodgson and Moore^. For less negative ax we 
use results for the correlation length directly from that 
paper. The fit is similar in appearance to Fig. [2] but 
we do not reproduce it here; it is slightly less convincing 
than the scaling in the 2D OCP plasma. 

We have found for the 2D OCP in the disc geometry 
that the radial density oscillations allow us to see directly 
the magnitude of the correlation length in the bulk fluid. 
The same behavior is found in the GL LLL model. It is 
now clear that those simulation studies which reported a 
first-order fluid-crystal transitionii^i^i^ involved systems 
whose linear dimensions at the reported transition were 
actually comparable with the correlation length. It would 
clearly be desirable to repeat these studies for larger sys- 
tems before accepting the conclusion that there is a first- 
order transition in these two-dimensional systems. 

Finally, we note that it should be possible to see the 
density oscillations we have reported in this paper at the 
edges of thin film type II superconducting samples. 
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